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I. INTRODUCTION 

A point-like particle diffusing in a two or three dimen- 
sional (2D, 3D) channel of varying cross section became 
an archetypal model describing transport through nano 
channels, pores or along fibers in biological systems, as 
well as passing of large molecules through membranes. 
Analytic studies of such models usually require further 
simplifications, namely the dimensional reduction to a 
purely one-dimensional (ID) system described by an ef- 
fective ID evolution equation, governing the linear (ID) 
density p{x,t), depending on time t and the longitudi- 
nal coordinate x. On the other hand, any simplification 
should retain all important features of the original full di- 
mensional model. After the dimensional reduction, they 
are reflected in the structure of the effective equation. 

The Fick- Jacobs (FJ) equation 



dtp{x, t) = D(^dxA{x)d, 



A{x) ' 



(1.1) 



can serve as the simplest example of such an effective ID 
equation for diffusion in a channel with reflecting walls; 
A{x) denotes the cross section area for 3D or the width 
for 2D channels at some x, and Dq is the diffusion con- 
stant. This equation maintains only the mass conserva- 
tion along the channel with varying A(x). Introducing a 
spatially dependent effective diffusion coefficient D(x) in 
the effective equation 0, 0| , 



dtp{x,t) ^ d^A{x)D{x)d, 



^ A{x) ■ 



(1.2) 



enables us also to respect boundary conditions (BC) and 
the local mass conservation at a point (x, y) of the full 
dimensional problem in the case of the stationary flow, 
i.e. when the net flux J{x, t) flowing through the channel 
is constant (y denotes the transverse coordinates). For 
nonstationary processes, one should replace the function 
D{x) by an operator D{x), containing also the spatial 



derivatives 9", n — 1, 2, but in practice, in most cases 
the asymptotic behavior of the processes is studied and 
the FJ equation extended only by the function D{x), Eq. 
()1.2|) . represents a significant improvement of the stan- 
dard FJ approximation (jl.ip . 

Of course, it is necessary to find a way how to fix the 
function D{x). Based on a fenomenological argumenta- 
tion, Reguera and Rubi Q suggested a function 



D{x) =Do[l + h'\x) 



(1.3) 



where rj — 1/3 or 1/2 for 2D or 3D channels with axial 
symmetry, respectively, and h{x) denotes half width or 
radius of the channel. Later Kalinay and Percus 
showed, that this function can be found as a series ex- 
pansion in a small parameter e, representing the ratio 
of the longitudinal and the transverse diffusion constant, 
e = Dq/ Dy. The anisotropy of the diffusion constant, im- 
posed artificially, causes separation of the modes quickly 
decaying in the transverse direction from much slower 
longitudinal ones and formally, it enables us to find a re- 
currence scheme, generating systematically higher order 
corrections to the F J equation (jl.ip , giving the expansion 
of the function D{x) in e in the stationary state. For 2D 
channels, bounded hy y = h{x) and the x axis, we get 



-h'^ + —h' 
3 45 



Dix) - i^o(l 

X [9ft'3 -t- hh'h" - 



(1.4) 



If h!'{x) and the higher derivatives are neglected, the re- 
sult is 



D{x)c±Dn[\ 



-h'' 



Do 



arctan y^h' 

7^' 



(1.5) 

a function differing from (|1.3|) by less than 1% for mod- 
erate slopes, \h'\ < 1, for an isotropic diffusion, e = 1. 
For 3D symmetric channels, the same treatment results 
in the formula (|1.3I) and rj — 1/2. All formulas exhib- 
ited good agreement with numerical tests [3] for \h'\ < 1; 
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FIG. 1: A sketch of the considered model: the channel is 
bounded by hard walls at y — h{x) and —h{x). Diffusing 
particles are biased by a constant gravitational force G 



steeper slopes require to take also the higher derivatives 
of h(x) into account |8|, or to avoid the expansion for 
specific geometries at all 13 1. 

Let us remark that introducing anisotropy of the dif- 
fusion constant is equivalent to rescaling of all transverse 
lengths (i.e. the coordinate y and the half width h{x) in 
the 2D channels) by y/e Q . This scaling together with a 
similar recurrence procedure were used to calculate cor- 
rections to the mean velocity and the dispersivity of the 
particles IJ, |15| within the macrotransport theory, as 
well as for re-derivation of D{x) in Eq. (|1.5p fl6'|. 
Recent studies 



17 



18| showed that the same strategy 
could be used also for the dimensional reduction of the 
Smoluchowski equation, 



dtp{x,y,t) 



DyVy ■ e-^(^'y)/'^-«^V,e^(-'y)/^-^]p(x,y,i),(1.6) 

i.e. for mapping of the diffusion in an external field 
U{x,y); p{x,y,t) is the 2D or 3D density, Vy denotes 
the gradient in the transverse directions and T is the tem- 
perature. The ratio e = Do/Dy, remains a good small 
parameter also for the biased diffusion, enabling us again 
to construct a recurrence procedure generating system- 
atically corrections to an equivalent of the FJ equation. 
This extension allows us to apply the dimensional re- 
duction to a much broader class of problems interesting 
in chemical physics. Considering a force parallel to the 
X axis of the channel, we showed [l7[ how the entropic 
potential is added to the real (energetic) one in the di- 
mensionally reduced dynamics, which can be useful for 
studying e.g. the Brownian pumps [l^. Mapping of Eq. 
p.6p for a potential depending on the transverse coordi- 
nates [l^ can be used to get the reduced dynamics in a 
channel with soft walls. 

In the present paper, we study diffusion in a 2D sym- 
metric channel, bounded by smooth functions y = h(x) 
and —h{x), with hard and reflective walls (Fig.l). The 
particles are biased in the transverse direction by a con- 



stant gravitational force G; the potential U(x, y) ~ Gy in 
Eq. (|1.6p . This model was investigated mainly in connec- 
tion with the stochastic resonance during the last years 
[20t - 24 1 . An important feature of this model is an inter- 



play between the gravitational force, holding particles in 
the potential wells in the wider parts of the channel, and 
the thermal motion, enabling the particles to diffuse into 
the neighboring compartments over the potential barri- 
ers formed by the narrowings of the channel. An oscillat- 
ing force applied along the channel can help this hopping 
very effectively at a specific resonance frequency, depend- 
ing on the force G, the temperature T and the geometry 
of the channel. 

We focus our attention on the competition between 
the gravitational potential and the "entropic" potential 
in diffusion through such channels. Our mapping proce- 
dure allows us to quantify their contributions to the net 
fiow of particles in an elegant way: in the form of an 
effective ID equation, involving both effects in its struc- 
ture. The previous analyses used the ID description, too, 
but governed by an equivalent of the FJ equation ()1.1|) . 
Comparison of this theory with the Brownian simula- 
tions 2^ 26 1 indicates that this approximation may be 
not satisfactory especially in the region of our interest, 
when both effects become comparable. 

In the following section, we present the rigorous di- 
mensional reduction of this model onto the longitudinal 
coordinate. We show how to merge the mapping of dif- 
fusion in a transverse field [3| with the presence of the 
refiecting hard walls. The result of our mapping is an 
equation of the type (jl.2|) in the limit of the stationary 
fiow, with D{x) expanded to the first few orders in e. 
In Section 111, we suggest and justify an interpolation 
formula for D{x) based on results of the mapping pro- 
cedure in the "linear approximation", when h"{x) and 
higher derivatives in the expansion of D(x) are neglected. 
Our formula is verified by an exactly solvable model, the 
stationary diffusion of particles in a linear cone. 



II. MAPPING PROCEDURE 

We follow the procedure developed for the diffusion 0, 
M and the biased diffusion Il8||, based on introducing 
a small parameter e into the 2D Smoluchowski equation 
()1.6|) . We define it as a parameter of anisotropy of the 
diffusion constant, Dy = Dq/ t, but it can be also imposed 
by scaling of the transverse lengths |16j , y — > ^y and 
the inverse scaling of the force G — )■ G/^. Anyway, we 
get the equation 



dtp[x, y, t) = dlp{x, y, t) + -dye-^ydye^yp{x, y, t) (2.1) 



from Eq. (|1.6p for the model of our interest; we rescaled 
time t by the diffusion constant Dq, DqI — > t, and g — 
G/keT. 
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The small parameter e enables us to carry out two 
important steps of the mapping procedure. First, we can 
find readily the equivalent of the F J equation in the limit 
e — > 0, and second, it becomes a parameter controlling 
the perturbation expansion of any quantity describing 
diffusion in the channel: the 2D density p{x, y, t) and 
the flux density j{x,y,t) [l^, ID density p{x,t), or any 
mean value, like the mean velocity or dispersivity [l^ . ll5| . 
representing a sequence of corrections to the zero-th order 
(FJ) solution. 

Before the mapping, we have to supplement boundary 
conditions (BC). The Smoluchowski equation (|2.ip rep- 
resents the mass conservation law, so the components of 
the flux density j are 

jx{x,y,t) = -dxp{x,y,t) , 

jy{x,y,t) = --^e-'^ydyeayp{x,y,t) . (2.2) 

No flux through the reflecting hard walls requires to have 
the vector j at the boundaries parallel to them, so we get 

e^y dye-ay p{x, t) ^ ±eh\x)dxp{x, y, t) (2.3) 

y=±h{x) 

at the upper and the lower boundary y = ±h(x). BC 
at the ends of the channel are arbitrary, they do not 
enter the mapping procedure in our formulation. We can 
consider the channel as inflnite. 

The mapping procedure reduces the 2D Smoluchowski 
equation (j2.ip governing the 2D density p{x, y, t) to some 
ID equation governing the ID density p{x,t), deflned as 



h{x) 

p{x,t) = / p{x,y,t)dy. 

-h(x) 



(2.4) 



Thus the first step of the mapping is integration of Eq. 
(12. ip over the cross section. Applying the definition (12. 4p 
on the left hand side, we arrive at 



dtpix, t) 



h{x) 



h(x) 



dlp{x,y,t)dy 



h{x) 



dx I dxp{x,y,t)dy 

'h{x) 



(2.5) 



after integrating by parts and using BC p.3p . 

Our goal is also to express the right hand side of Eq. 
()2.5p in terms of p{x, t) instead of p{x, y, t). This task is 
easy to complete in the limit e — > 0. For an infinitesimally 
small e, the transverse diffusion constant Dy becomes al- 
most infinite and the 2D density po is equilibrated in the 
transverse direction almost immediately after any change 
in the x direction. So we can write 



Pa{x,y,t) 



1 



Aix 



■e yyp{x,y), 



(2.6) 



where A{x) provides normalization of po- If substituted 
in the condition (j2.4l) . we have to obtain an identity. 
Hence 

A{x) = I ^ ^ e-yydy = - sinh [.g/i(a;)]. (2.7) 



If we use the formula (|2.6p for the 2D density in Eq. (12.51) , 
we find 



dtp{x,t) 



h{x) 



-gy 



dy 



-h(x) 



dxA{x)dx 



p{x,t) 
A{x) ' 



(2.8) 



which is (an equivalent of) the FJ equation (|l.ip . Let us 
stress that A{x) is not the width of the channel here, but 
the integral (|2.7p . On the other hand, for 5 — >■ 0, A{x) 
becomes 2h{x). 

For e > 0, the transverse diffusion constant Dy is fi- 
nite and the local equilibrium in the y direction is dis- 
turbed by the flux flowing along the curved boundaries 
at a given x. So the formula ()2.6p cannot be used for the 
2D density, po{x,y,t) does not satisfy the Smoluchowski 
equation (12. ip . The small parameter e enables us to look 
for deviations of the real 2D density p{x,y,t) from the 
equilibrated density po l|2.6p in the form of a sequence of 
corrections, as an expansion in powers of e. 

Instead of expanding some speciflc density p{x,y,t) in 
e, we prefer a more general method. We look for the 
expansion of a wide class of solutions of the 2D problem; 
we expand all p{x^ y, t), which can be generated from any 
ID solution p{x, t) of the searched ID equation by the 
backward mapping onto the space of solutions of the 2D 
problem. Formally, these 2D densities can be expressed 
by the formula 



p{x,y,t) = e <^yu}{x,y,dx) 



A{x) ■ 



(2.9) 



where uj represents an operator of the backward mapping. 
It acts on a wide class of p{x,t), so the dependence of 
p{x,t) on e is not important, we take it as an unfixed 
function. Instead, we look for expansion of the operator 



UJ m e, so 



p{x,y,t)=e yy^e'''u}n{x,y,dx) 



71=0 



A{x) ■ 



(2.10) 



We know already the zero-th order term, uj^ — 1, which 
gives the equilibrated solutions po (|2.6I) in the limit e — ?> 0. 

Supposing p(x, y, t) of the form (|2.10p . we can formally 
complete the construction of the ID evolution equation. 
Applying this expression in Eq. (j2.5l) . we find 



dtp{x,t) = dx I dye-s^9, Ve"c2;„(x,y,a,)i^ 



4 



A{x) ■ 



(2.11) 



Here we already applied Wq = 1 and introduced an oper- 
ator Z, correcting the FJ equation (|2.8I) . It can also be 
expanded in e, 



Zk{x,dx)dx = - 



fc=i 
1 



h(x) 



h(x) 



dye-syd^LOk{x,y,d^). (2.12) 



Finally, the 2D density ()2.10p has to be a solution of 
the original Smoluchowski equation p.ip . 



3=0 



[e-^y{dt~dl)--^dye~^ydy]uj,{x,y,d,)P^ = 



Aix) 



(2.13) 

which generates a recurrence relation fixing the operators 
ujn- Because we suppose that these operators act only 
on the spatial coordinates, the time derivative commutes 
with them and for dtp{x,t), we use the equation (j2.11l) . 
Collecting the terms at the same powers of e, we find 



dye yydyLUn+i{x,y,dx) ^ - 

n ^ 

k=0 ^ ' 



-gy 



dlujn{x,y,dx) 



d.,A{x)Zk{x,d^)d^ ; (2.14) 



we take Zo{x,dx) = —1 in this formula. After double 
integration, we obtain uJn+ij giving the n+l-st correction 
to the 2D density, if applied on some ID solution p{x, t). 
Two integration constants have to be fixed (they are also 
operators, but independent of y). The first one provides 
satisfaction of the BC (|2.3I) . Using the formula (|2.10p in 
Eq. ()2.3p and comparing the terms at the same powers 
of e, we get the condition 



dyU)n+i{x,y,dx) = ±h'{x)dxUJnix,y,dx) 



y=±h{x) 



(2.15) 

If the integration constant is fixed at one boundary, the 
BC at the opposite boundary is automatically satisfied. 
The second integration constant helps to satisfy the nor- 
malization condition; applying the formula ()2.10p in the 
definition (|2.4[) has to give identity in any order of e, 
hence 



h{x) 



dye yyLdn{x,y,dx) = 



(2.16) 



-h(x) 



for n > 0. By this condition, we keep the same number of 
particles in both, ID and 2D descriptions, independently 
of how the densities p ov p are normalized. 

The recurrence procedure starts from the FJ approx- 
imation. Wo = 1 a-iid Zq = —1. Having expressed uJn+i 
for some n, the next order correction operator Z„+i is 



calculated according to the formula (|2.12p . Calculation 
of the first order correction and other details are given in 
the Appendix A. We show here only the results. 



h' 

UJl = — 

9 



+ (1 - 5y) cosh gh 



sinh gh 



-ghil 



sinh gh 



dx 



and the corresponding 
h'^ 



Zi 



sinh gh 



1 + cosh gh — 2ghcothgh 



(2.17) 



(2.18) 



One can check that in the limit 5 — >■ 0, we obtain Zi — 
h^/S, known for the diffusion alone 0, Q- The higher 
order operators Z„ starting from n — 2 also contain the 
spatial derivatives dx, what makes the equation (12. lip 
too difficult for direct use in practice. Alike the diffusion 
alone Q , this equation can be simplified by replacing the 
correction operator 1 — eZ by the function D(x) in the 
limit of the stationary state, when the net flux changes 
very slowly. 

In that case, Eq. (12. lip is replaced by an equation of 
the form p. 21) . where A{x) is given by the formula (|2.7p 
and D{x) has to be fixed. Thus we have two different 
expressions for the net flux. 



J{x,t) = -A{x) l-eZ{x,dx) 



ft 



^ A(x) 



and 



J{x,t) = -A{x)D{x)d, 



Aix) ■ 



(2.19) 



(2.20) 



coming from Eqs. (12.111) and (|1.2p . respectively, as both 
equations represent the ID mass conservation law. In 
the stationary state, J{x,t) = J is constant in time and 
space and dx[p{x,t)/A{x)]/ J = —l/A{x)D{x) depends 
only on geometry and the parameters of the model for 
any stationary solution p{x,t) = p{x). Then the formula 
(|2.19p describes the same flux J only if 



1 



D{x) 



A{x) l-eZ{x,dx) 



1 



A{x) ' 



(2.21) 



which fixes the effective diffusion coefficient D{x) unam- 
biguously for Z obtained from the mapping procedure. 
If the expansion of Z in e (|2.12p is used in Eq. (|2.2ip . 
the result is an e-expansion of D{x), 



D{x) = 1- 



eh'- 



+- 



sinh gh 
e^h'^ 



sinh gh 



1 + cosh gh — 2gh coth gh 
gh 



sinh gh cosh gh ^ smh{2gh) 



X (17 sinh^ gh + 36) + {gh f (7 sinh"* gh 
+40sinh2 5/i + 36) + O(e^) + 0(/i"); (2.22) 

the second and higher derivatives of h{x) are already ne- 
glected in this formula. 
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FIG. 3: Linear cone with a constant transverse force G. 



FIG. 2: The coefBcient D(x) depending on the local slope 
h'{x) and the values of gh{x) = 0,0.5,1,2 and infinity at 
some point x. The thick lines depict the limits (/ — >■ and oo. 
The dashed lines correspond to the interpolation formula (|3.2p 
with the exponent (|3.3p . The adjacent thin full lines describe 
the truncated expansion (|2.22p up to the 3-rd order (the lower 
lines) and the 4-th order (the upper lines). The dots depict 
the data gained from the exactly solvable model, a linear cone 
with h'{x) = tan(7r/10) ~ 0.325 and tan(7r/6) = 1/^3. 



III. INTERPOLATION FORMULA FOR D(x) 

Even in the " linear approximation" , which neglects all 
but the first derivative of h{x) in the expansion of D{x), 
the resulting formula (|2.22p is much more complicated 
than the similar one valid for the diffusion alone, Eq. 
(ll.Sp . Also it is difficult to sum directly the series in 
eh'^ up to infinity, to find a formula for D{x) in a closed 
form, usable in practice. In this section, we propose an ad 
hoc formula, justify its validity and test it on an exactly 
solvable model. 

The expansion (I2.22p simplifies considerably in two 
limits: for g — > and gh{x) — > oo. The first case 
corresponds to the unbiased diffusion, the coefficients 
at (e/i'^)" approach (— l)"/(2n -I- 1) and the series is 
summable, giving the result (|1.5p . In a strong gravita- 
tional field, the limit of the coefficients is (—1)", so 



D{x) 



1 - eh 



12 



2il4 



1 



1 + e/i'2 ' 



(3.1) 



the proof is given in the Appendix A. 

Finally, we recall that the formula (11.31) differs only 
slightly from the exact result, Eq. (|1.5p . at moderate 
slopes of the walls, |/i'(2;)| < 1. Then it seems reasonable 
to suggest this formula also for the region of intermediate 
g, but with the exponent —77 depending on gh{x), 



D{x) c^Dq[1 + eh'^{x)]-'i^sh{x}] 



(3.2) 



For the choice 



1 



ri[gh{x)] = - o 

smli gn 



I + cosh^ gh — 2ghcoth gh , (3.3) 



we recover correctly the first order term of D(x) in Eq. 
(12221) ■ Then in strong fields, the formula ([5^ ap- 
proaches the exact limit p.ip and for g — J- 0, we get 
the function of Reguera and Rubf, Eq. (|1.3p . 

We compare first our interpolation formula with the 
true expansion of D{x) (|2.22p calculated up to the 4- 
th order in e. The plots of D{x) versus slope of the 
boundaries h'{x) are depicted in Fig. 2. The thick lines 
describe the limits, g — >■ and 00. The dashed lines 
plot the formula p.2p for three intermediate values of 
g = 0.5, 1 and 2. These data are compared with the 
truncated series (|2.22l) . the adjacent thin lines include 
the corrections up to the 3-rd order (the lower lines), 
and the 4-th order (the upper lines). In the region of 
fast convergence of the series (I2.22p . where the lines of 
the 3-rd and the 4-th order formulas almost coincide, the 
difference between the true and the interpolated values 
is comparable with the difference between the formulas 
([T3)) and do]) (the thick dashed line). 

Unfortunately, the radius of convergence of this series 
is finite and decreasing with growing g. So we test our 
interpolation formula on an exactly solvable model. 

Tests of such theories are often based on calculation 
of the net flux J flowing through an exactly solvable 
structure. The flux calculated from the exact 2D den- 
sity p{x,y), 

/h(x) r^{^) 
jxix,y)dy = - dxp{x,y)dy (3.4) 

'h(x) J-h{x) 

in the stationary regime, is compared with the corre- 
sponding flux according to Eq. (|2.20p with D{x) derived 
within the tested theory. We modify this method: for a 
given exact solution p{x, y), we calculate the flux J (13.41) . 
the ID density p{x) (|2.4p and the corresponding D{x), 



D{x) 



J 



dx 



p{x) 



A{x) \ A{x) 



(3.5) 



from Eq. (j2.20p . which is compared with D{x) coming 
from the theory, Eq. p.2p or ()2.22p in our case. 

Our exactly solvable model is a stationary flow through 
a linear cone, bounded hy y — zLax, see Fig. 3. The 
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FIG. 4: Contour plot of the 2D density p{x, y) (the left panel) 
and p{x,y) — e^^p{x,y) (the right panel) according to Eq. 
p.6|) in a channel bounded hy y = ±a;/\/3 and g — 2. 



particles are emitted from a point-like source at the origin 
of the coordinate system and collected at an absorbing 
boundary placed far from the positions x of our interest. 
Let us notice that the expansion (|2.22p of D{x) summed 
up to infinity describes our model exactly; h'{x) = a is 
constant and its derivatives are zero. 

In the Appendix B, we show that the 2D density ex- 
pressed in the form of an integral in the complex plane 



r'tTT /2-\-oo 

-gy/2 I p-ig^/x'^+y^/'^) cosh(2-j7r/2) 



X f{z + + f{z — icj) ~ in) dz + c.c, (3.6) 



(j) — arctany/x and f{w) ~ coth(mw/2) tanh(i(;/2), m — 
3,5,7, represents a stationary solution of the Smolu- 
chowski equation (12.11) with BC (|2.3p . e — I, Dq — I and 
h(x) = ax for specific values of the slope a — tan(/)o; 
00 — 7r/2TO = 7r/6, 7r/10, ... The integration from to 
i7r/2 -t- oo (and to — i7r/2 -t- oo in the complex conju- 
gated expression) is carried out along any path avoiding 
the poles of the integrand on the imaginary axis from 
the right side, see Fig. 7 in the Appendix B. The con- 
tour plots of the density p (|3.6p and the corresponding 
p{x,y) = e^^p{x,y) are shown in Fig. 4. According to 
Eq. (|2.2p . the gradient of p{x,y) is proportional to the 
flux density, so one can check visually on the right panel 
that the no flux BC are satisfled on both boundaries. 

For testing the interpolation formula (13. 2p . we use the 
channels with a — tan(7r/6) = l/VS and tan(7r/10) ~ 
0.325. The values of p{x) and J in Eq. p.6p were inte- 
grated numerically; the calculation of J serves as a test 
of the numerical method, since J does not depend on x. 
A choice of the force g is not important; it scales the 
length unit in both directions, as can be seen from Eq. 
(|3.6p . Finally, we express a and x by using h{x) and 
h'{x), a = h'{x) and x = h{x)/h'{x), valid for the linear 
cone, to place the results in the plot of D{x) depending 
on h'{x) and gh{x). 

The data for gh{x) — 0.5, 1 and 2 are depicted as 



1 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 



2 

gh(x) 



FIG. 5: The exponent t) plotted versus gh{x) according to Eq. 
(13. 3p (the line) and gained by fitting the interpolation formula 
p.2p to D{x) calculated for the linear cone, 4>o — tt/G (the 
larger dots) and tt/IO (the smaller upper dots) 



dots in Fig. 2. The interpolation formula describes 
the coefficient D{x) satisfactorily for small slopes of the 
boundaries h' and close to the limits gh{x) — and oo. 
For larger h' > 0.5 in an intermediate region, roughly 
1 < gh{x) < 5, the deviations are more notable. For 
practical purposes, one could try to find an interpolation 
formula for rj fitting better the exact data obtained for 
the linear cones. The exponents of Eq. p.2p fitted to 
the exact values of D{x) for the cones with a = tan(7r/6) 
(the larger dots) and tan(7r/10) (the smaller dots) are 
depicted in Fig. 5 and compared with the function 



IV. CONCLUSION 

The main aim of this paper was to arrive at an effective 
ID description of diffusion in a 2D symmetric channel of 
varying width 2h{x). The diffusing particles are biased 
by a constant gravitational force G acting in the direction 
perpendicular to the axis of the channel. 

Our effective equation of the type ()1.2p , governing evo- 
lution of the ID density p(x, t) in the channel, goes be- 
yond the Fick- Jacobs approximation, considering only in- 
stant equilibration of the 2D density in the transverse di- 
rection, which was used in the studies based on this model 
2ol 25 , 2^ till now. The effects of slower transverse re- 
laxation are included in the effective diffusion coefficient 
D{x). We calculate this function within a recurrence pro- 
cedure mapping rigorously the 2D problem onto 
the longitudinal coordinate x in the limit of the station- 
ary flow, i.e. when the net flux changes very slowly with 
respect to the relaxation in the transverse direction. 

The result is an expansion of D{x) (|2.22p in a parame- 
ter e expressing the ratio of the diffusion constant in the 
longitudinal and the transverse directions, e = Do/Dy, 
introduced artificially in the Smoluchowski equation (|2.ip 
and set to 1 at the end. Adding the transverse force 



7 



makes the result much more comphcated, if compared 
with D{x) (11.51) for the diffusion alone. It is difficult to 
obtain a simple formula usable in practice by direct sum- 
ming of the expansion in eh''^(x) up to infinity. So we 
suggest to use the interpolation formula p.2p introduced 
before by Reguera and Rubi [3|] for diffusion. We showed 
that the biasing force effectively changes the exponent 77, 
depending on gh{x), g = G/ksT. It increases from 1/3 
for negligible G up to 1 in strong fields. This dependence 
can be approximated by the function (13.31) ; then the first 
order correction of the exact D{x) (|2.18p is recovered. 

The interpolation formula is compared with the trun- 
cated exact expansion (I2.22p up to the 4-th order and also 
checked by the model of biased diffusion in a linear cone, 
which is exactly solvable. The agreement is satisfactory, 
but the differences increase especially at steeper slopes 
of the boundaries, h'{x) > 0.5, and intermediate values 
of gh{x) (roughly units). For practical purposes, one can 
try to find some other (ad hoc) function for the expo- 
nent r]{gh{x)) to fit better the exact D{x) in this region. 
Let us recall that the "linear approximation" does not 
work well for h' > 1 [3], because higher derivatives play 
an important role there, too. Then taking more com- 
plicated interpolation formulas depending only on h'{x) 
may have only a small effect on the further improvement 
of the results. 

The calculation of the expansion of D{x) ()2.22p pre- 
sented in the Section II also demonstrates how the map- 
ping procedure [1, @| can be applied to diffusion bounded 
in a channel with hard walls and biased by a transverse 
force. Other possible extensions are straightforward: we 
can add also a force acting along the channel, or to go to 
3D channels. The problem is growing complexity of the 
expansions in e and necessity of summation of at least 
some group of their terms up to infinity; the expansions 
are converging only in a restricted region of parameters, 
as seen in Fig. 2. An effective way of searching for D{x) is 
combination of the mapping recurrence scheme with fit- 
ting the results of exactly solvable models, as suggested 
in [8| and applied also in this work. 
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APPENDIX A: DETAILS OF MAPPING 

We demonstrate here the mapping procedure on cal- 
culation of the first order correction and then we prove 
the formula p.l|) in the limit of large g. 

Starting from the zero-th order, we take n = in the 
recurrence relation (|2.14p . ujq ^ 1 and Zq = —1. We 



obtain 



a, 



dyCoi = e-^y -rd^Ad, -dt\= e-sy—d, 



A 



.A' 



A 



(Al) 



After applying Eq. (|2.7p and the first integration, 

dyUJi = e^y j dye-3ygh\x)cot\Ygh{x)d^ 

-h'{x) coihgh{x)e~yyd.:, + Ci] , (A2) 



we fix the first integration constant Ci from the BC (I2.15P 
at the upper boundary, dyUii ~ h'{x)dx at y — h(x), 

Ci = /i'(x)e-s''(=") [1 -I- coth gh{x)] d^. (A3) 

Notice that Ci is an operator, but independent of y. Also 
one can check that the relation 



dyUJl = — 



h'{x) 



s'mh gh{x) 



[e^y — coshgh{x)] dx 



(A4) 



satisfies the BC (|2.15l) at the lower boundary, y — —h{x), 
too. The next step is integration of Eq. (jA4p . 



h'{x) 



sinh gh{x) 



-esy -y cosh gh{x) + Cq, (A5) 
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and fixing the second integration constant Cq (again an 
operator) from the normalization condition ()2.16p . After 
some algebra, we arrive at the formula (j2.17p . Finally, 
applying Eq. (|2.12p to the resultant uji gives the first 
order correction operator Zi (|2.18l) . 

In the limit g ^ 00, we keep only the leading terms of 
any expression during the calculation; the other terms, 
proportional to powers of e~^^^^\ are negligible. For 
A{x) ~ (l/g)e^'''(^\ the initial equation ([M]) of the re- 
currence scheme becomes 



dye-yydyUJi = gh'{x)e~yydx ■ 



(A6) 



After the first integration and fixing Ci at the upper 
boundary, we get 

dyibi = h'{x) (2eS(^-^(-)) ~l)dx; (A7) 

BC are satisfied at the lower boundary, too, because the 
term ^ e~2g/t(a;) -g negligible. The second integration and 
fitting the normalization condition gives 



LUl = h'(x) 



9 



9 



dx 



(A8) 



up to the terms ^ e"^''^^^ and smaller. Finally, in the 
integration of Zi according to Eq. (|2.12p . only one term 
remains. 



Zi = ge 



h{x) 



h{x) 



dye-yyh'^{x) 



(A9) 



8 



nonvanishing in the limit of large g. Notice also that 
the exponential term in Eqs. (jA7[) and (|A8|) does not 
contribute to Co and Zi] we can neglect it. This simpli- 
fication corresponds to fixing Ci at the lower boundary, 
then dyUJi — —h'{x)dx- The BC at the upper boundary 
is not satisfied, but on the other hand, there are no par- 
ticles there for large g. We treat the upper boundary like 
it was in infinity. 

Now we can prove the formula p.f p . First we simplify 
the equation (|2.2fp . 



1 



D{x) 



A{x) 



1 + eZ 



1 



Aix) 



1 + + 



1 2 



eA{x)Z 



Aix) 



,,(AfO) 



hence 

D{x) ~ l-eA{x)Z 



OO 



A{x)' 
(All) 

the difference depends on derivatives higher than h'{x) 
and they are neglected in our "linear" approximation. 

The terms of the series in Eq. (|Alip can be ex- 
pressed directly by applying the operators c2;„ on a func- 
tion f{x) = J dx/A{x) ~ / ge-a^^'^Hx. Then 



A{x)Z„ 



A{x) 



h{x) 



h(x) 



dye-ay dxLo^f{x) (A12) 



from the relation (|2.12p . The functions c2;„/(x) are de- 
rived by the same recurrence procedure, as it was demon- 
strated on the operator Coi above. If we retain only the 
leading terms in the limit f/ — >■ (X) in our expressions and 
neglect h"{x) and its derivatives, we arrive at 

w^fix) - -h'{x)e-a'''^-^giy + hix))~l], 
^2f{x) = -—e-<^^^[g\y + hf-Ag{y + h) + 2\, 

ojsfix) = -l!^e-<^'^[g^y + hf^9g\y + hf 



+18g{y + h)-6], 



k=0 



kl 



[giy + ht. 

(A13) 



We can check normalization (|2.16p of these formulas, 



dye 



E 



(-1) 



n — k 



E 

fc=0 



(-1) 



n — k 



kl 



«]' 



9 kl 



01 



e''^z^dz = 0, (A14) 



after substituting z = g{y + h) and replacing the upper 
boundary 2gh by infinity; we omitted writing explicit de- 
pendence of h(x) on X. The coefficients of the expansion 



(jAlip are integrated in a similar way; after completing 
the X derivative of the formulas (jA13p in Eq. (jA12p and 
using the normalization (IA14p . only the term 

A{x)Zn 



A{x) 



= /i'2« = (A15) 

fc=i 



remains, the result we wanted to prove. 

Finally, one can check by direct calculation that the 
formulas (jA13p satisfy the recurrence relation (|2.14l) and 
the BC (I2.15P at the lower boundary. The operators w„ 
are replaced here by the functions c2;„/(x), dxf{x) — 
g^-gH^)^ The terms depending on Z^. disappear from 
Eq. (EUl), since A{x)Zk{l/A(x)) = (_i)("+i)/i'2n ac- 
cording to Eq. (|A15p . and its derivative depends on h'\ 
which is neglected. 



APPENDIX B: EXACT SOLUTION 

We present here the stationary solution of the Smolu- 
chowski equation (|2.ip for the biased diffusion in a linear 
cone. 

For a point-like source of particles placed at the origin 
of the coordinate system, the stationary equation (|2.ip . 

= dlp{x, y) + dye-^ydye^ypix, y) (Bl) 

becomes separable after substitution 



-gyf^ 



u{x,y) 



(B2) 



and converting to the polar coordinates, x — r cos 0, y 
r sin (h. 



u(r,(/))=0. (B3) 



The particular solutions are u{r,(f)) = i?„(gr/2)e"'"^; i?„ 
stands for the Bessel /„ or the Bessel functions. In an 
unbounded plane, the linear combinations describing the 
biased diffusion cannot contain the /„ functions, since 
p{x, y) would diverge for y — > ~oo. The solution u(x, y) 
is then composed from the Bessel Kn functions. The 
"ground-state" solution ua{r,(j>) = Ko{gr/2) carries the 
flux along the force. Other particular solutions u„(r, 0) = 
Kn{gr/2) sm{n(f)) are necessary for fitting p(x,y) = at 
an absorbing boundary, if it is considered, but they do 
not contribute to the net flux; one can check that 



u„ix,y)eyy/^ 



dx = —JS„ 



(B4) 
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FIG. 6: The linear cone in Fig. 3 rotated by 7r/2 



at any fixed y < (below the source). We need only 
the ground state uq{x, y) for the calculation of D{x), Eq. 

Our problem is to find such a solution for diffusion in 
a linear cone, i.e. satisfying the BC (12.31) at ?/ = ±ax. In 
the polar coordinates, the BC become 



£?0u(r, 4>) 



gr 



cos{(f>)u{r, (p) 



i>=±4 



(B5) 



for any r > and h'(x) ^ a — tan^Q. 

This task is related to the calculation of the 2D sta- 
tionary density of particles dragged out of the cone by 
a constant force along the x axis [l3|. If we rotate our 
channel in Fig. 3 by 7r/2, we get the relevant picture. 
Fig. 6. In comparison to the previous problem, the 
particles diffuse in a different sector; instead of the an- 
gle ip = <j) + 7r/2 e (0,7r/2 — ^o), they are confined in 
tp € (7r/2 — (j)o,TT/2 + (j)o). For certain values of (po, we 
can extend the known solutions in the sector adjacent to 
the X axis [17] to the sector of our interest. 

We recall briefly the stationary solution of the Smolu- 
chowski equation in the sector ■0 g (0, ipo). After rotation 
of the coordinate system, (p is simply replaced by tp in Eq. 
(jB3| and the rotated BC (jB5|) . 



(B6) 



gv 

d^u{r,ip) = -—sin{ip)u{r,%p) 



has to be satisfied dX ip — ±V'o- To express the solutions 
u{r,ip) here, we are inspired by the integral representa- 
tion of the Bessel functions K^, 27|, 
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-''""'^^ cosh I't dt . (B7) 
One can check by direct calculation, that the integral 



ii(r, "ip) — I e 
'o 



' {gr 12) cosh t 



f{t + tij) + f{t-i^P) 



dt , 
(B8) 

is a solution of the equation (|B3I) (with (p replaced by 
Ip) for any even analytic function f(z) — f{—z) of the 
complex variable z having no pole along the integration 
path. Notice also that the function (jB8|) has expected 



symmetry u{r,ip) — u{r^—ip), given by the direction of 
the force along the x axis. 

The function f{z) is fixed from the BC (jB6p ai 'tp ^ 
±^Po- Applying the formula (IBSp in Eq. (|B6p and inte- 
grating by parts, we obtain the condition 

[f{t + t^o)- fit-ii'o)] sinht- 

= t[f{t + i^/jo) + f{t~ ^V'o)] sin V-o (B9) 

valid for any t > 0. If we write f{z) — ^(z) tanh (z/2), 
g{z) has to satisfy g{t + iipo) — g{t — iipo). Then the 
"ground state" solution uo{r,ip) is generated by go{z) — 
cot]i(Trz/2ip()) and the other particular solutions M„(r, -0) 
come from gn{z) = sinh(n7r2/'0o)- Again, only uo(^j "0) is 
connected with the ID stationary flux flowing along the x 
axis, and Wn^g are modes projected out by the mapping 
procedure 17], which are not necessary for calculation of 
D{x). 

To get the solution uq for the cone with the transverse 
field, i.e. for the sector ip S (7r/2 — 0o, 7''/2 + (^o), we 
have to find the function f{z) such that the BC (|B6p are 
satisfied at ip — 7r/2 ± 0o- The same treatment leads 
to a condition similar to Eq. (jB9P : if we write f{z) = 
g{z) tanh(2/2), then g{t -t- iip) = g(t — iip) is required at 
both boundaries, ip = 7r/2 ± (pQ, and any t > 0. 

For specific angles (po, we can adopt the function 



go{z) = coth(mz/2) = 



1 



1 



(BIO) 



It satisfies the required condition not only at ipo = 7r/m, 
used in the sector adjacent to the x axis, but also at any 
its integer multiple. To get the "ground state", we need 
to adjust 7r/2 ± (po to be succeeding integer multiples of 
tt/to; the imaginary part of m{t + iip) has to change by 
in if Ip increases from tt/2 — cpQ up to Tr/2 + (pQ. These 
requirements are met for odd numbers m > 3. For the 
corresponding angles cpQ = 7r/2m = 7r/6,7r/10, .., the for- 
mula (jBlOp becomes the function go{z) generating the 
"ground state" ito(r, ip) also in the sector of our interest. 

Still, there is a problem at the boundary whose angle 
Ip is an even multiple of ipo] the function go{z) and also 
the corresponding f{z) have a pole a.t t — 0. We solve it 
by changing the integration path in the complex plane. 

First we return back to the unrotated coordinate sys- 
tem and the angle (p. The variable t in the integral (IBSP , 



u{r,(p) = I e 

'0 



-(3r/2)cosht J{t + icP + iTT/2) 

+ f{t-i(p~i7T/2)]dt , (Bll) 



can be substituted by t = z ± ?7r/2 and the path is then 
shifted in the complex plane by =pi7r/2 correspondingly. 
In the final formula, rewritten in a symmetric way. 



2u(r, cP) 



i-K /2-\-oo 



-{gr /2) cosh.{z— i-K 1 2) 
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FIG. 7: Integration paths in Eq. (|B12|) . The paths crossing 
poles on the imaginary axis (the dashed hnes) are replaced by 
the full lines, avoiding the poles from the right. 



+/(z — i(f> — itt) dz 

-i7r/2+oo 

-(gr/2) cosh(^+i7r/2) 



dz, 



f{z^i4>) 

(B12) 



we change the lower limits ±j7r/2 by zero and the integra- 
tion path in both integrals avoids the poles at z = ±in/2 
from the right, see Fig. 7. 

Notice that the function (|B12p holds the symmetry 
u(r,(j)) = u{r,—(j) — tt), determined now by the y di- 
rection of the force (Fig. 3). Direct calculation shows 
that also this u(r, 0) with changed lower limits of inte- 
gration solves the equation (|B3|) for any f{z) = f{—z), 
which has no pole along the integration path. Substitut- 
ing Eq. (jB12p in the BC (jB5|) and integrating by parts 
results in the conditions g{z ± = g{z =F =F m) at 
(j) = ±00 and any z on the integration path; we rewrote 
again f{z) = g{z) tanh(z/2). It is easy to verify that the 
function (jBlOl) satisfies these conditions for 00 — 7r/2m, 
m = 3, 5, so taking 



f{z) = coth(mz/2) tanh(z/2) 



(B13) 



in Eq. (|B12[) . we obtain the "ground state" solution 
Uf){r, 0) in the linear cone with a transverse force for these 
specific slopes a = tan0o- 

In the resultant stationary density, two integration 
constants can be added, 



p{x,y) = Cie-f^/2wo(x,2/) -f Coe-3^ 



(B14) 



Ci controls the stationary net flux connected with the 
density p{x, y) and Cq sets the BC p{x, j/) = at a dis- 
tant boundary absorbing the particles. None of them 
influences the calculation of D{x). The contribution of 
the term proportional to Co to p{x) / A{x) is constant and 



so it gives zero in the formula p.Sp . The net flux J, as 
well as dx[p{x) / A{x)], are proportional to Ci and so it is 
canceled in the resultant D{x). 
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